subroutine share_space_x
use var
implicit none
	integer i,j

	real(8) b0,b1,b2,bm1!,fpp(-bl:mx+bl),fmm(-bl:mx+bl),fp,fm
	real(8) q0,q1,q2,q3
    real(8) eps
	real(8) beta0,beta1,beta2,beta,beta3
	real(8) w0,w1,w2,w,w3
	real(8) Cui
	real(8) s ! for MUSCL
	real(8) theta,is2 !for weno3is+up3 ACFD11
	real(8) is51,is52 !for case 47
	real(8) is1,theta1,theta2 !for hybrid central4+weno3is
	real(8) f10,f11,f12,f20,f21,f22,g1,g2 !for WENO5
	real(8) tao !for WENO-Z
	real(8) qa !for weno5 paper
	real(8) xi !for WENO-NS
	real(8) c400,c401,c410,c411,c50,c51 !for multi-weno shen
	real(8) tao40,tao41,h0,h1
	real(8) alpha0,alpha1,alpha2,alpha3
    real(8) fs1(0:2),fs2(0:2)
    integer qw
    real(8) Ma(1+Klx:gridsx+krx,1+Kly:gridsy+kry),share_f(1+Klx:gridsx+krx,1+Kly:gridsy+kry)

	select case (c_case)
    case(9,13,16) !
        eps=1E-40
        share_f(1+Klx:gridsx+krx,1+Kly:gridsy+kry)=d(1+Klx:gridsx+krx,1+Kly:gridsy+kry)*p(1+Klx:gridsx+krx,1+Kly:gridsy+kry)*Qmean(1+Klx:gridsx+krx,1+Kly:gridsy+kry,4)      
    case(10,14) !
        eps=1E-40
        Ma(1+Klx:gridsx+krx,1+Kly:gridsy+kry) = (u(1+Klx:gridsx+krx,1+Kly:gridsy+kry)**2+v(1+Klx:gridsx+krx,1+Kly:gridsy+kry)**2)/c(1+Klx:gridsx+krx,1+Kly:gridsy+kry)**2
        share_f(1+Klx:gridsx+krx,1+Kly:gridsy+kry)=d(1+Klx:gridsx+krx,1+Kly:gridsy+kry)*p(1+Klx:gridsx+krx,1+Kly:gridsy+kry)/(Ma(1+Klx:gridsx+krx,1+Kly:gridsy+kry)+1.d0) 
    case(15) !
        eps=1E-40
        !share_f(1+Klx:gridsx+krx,1+Kly:gridsy+kry)=d(1+Klx:gridsx+krx,1+Kly:gridsy+kry)*p(1+Klx:gridsx+krx,1+Kly:gridsy+kry)*Qmean(1+Klx:gridsx+krx,1+Kly:gridsy+kry,4)  
        share_f(1+Klx:gridsx+krx,1+Kly:gridsy+kry)=d(1+Klx:gridsx+krx,1+Kly:gridsy+kry)+p(1+Klx:gridsx+krx,1+Kly:gridsy+kry)+Qmean(1+Klx:gridsx+krx,1+Kly:gridsy+kry,4)  
    case default
        stop'share_space_x: zero pivot at first row'
    end select
    select case (c_case)
    case(9,10,15) !三阶weno-z，upwind
        do j=1,gridsy
            do i=0,gridsx
                b0=share_f(i,j)-share_f(i-1,j)
                b1=share_f(i+1,j)-share_f(i,j)
                beta0=b0*b0+eps
                beta1=b1*b1+eps
                tao=dabs(beta0-beta1)
                beta0=(beta0+tao+eps)/(beta0+eps)
                beta1=(beta1+tao+eps)/(beta1+eps)
                w=1.d0/(1.d0+2.d0*beta1/beta0)
                !
                b0=d(i,j)-d(i-1,j)
                b1=d(i+1,j)-d(i,j)
                Qx_m(i,j,1)=d(i,j)+(w*(b0-b1)+b1)/2.d0
                !
                b0=u(i,j)-u(i-1,j)
                b1=u(i+1,j)-u(i,j)
                Qx_m(i,j,2)=u(i,j)+(w*(b0-b1)+b1)/2.d0
                !
                b0=v(i,j)-v(i-1,j)
                b1=v(i+1,j)-v(i,j)
                Qx_m(i,j,3)=v(i,j)+(w*(b0-b1)+b1)/2.d0
                !
                b0=p(i,j)-p(i-1,j)
                b1=p(i+1,j)-p(i,j)
                Qx_m(i,j,4)=p(i,j)+(w*(b0-b1)+b1)/2.d0
    
                b0=share_f(i+1,j)-share_f(i+2,j)
                b1=share_f(i,j)-share_f(i+1,j)
                beta0=b0*b0+eps
                beta1=b1*b1+eps
                tao=dabs(beta0-beta1)
                beta0=(beta0+tao+eps)/(beta0+eps)
                beta1=(beta1+tao+eps)/(beta1+eps)
                w=1.d0/(1.d0+2.d0*beta1/beta0)
                !
                b0=d(i+1,j)-d(i+2,j)
                b1=d(i,j)-d(i+1,j)
                Qx_p(i,j,1)=d(i+1,j)+(w*(b0-b1)+b1)/2.d0
                !
                b0=u(i+1,j)-u(i+2,j)
                b1=u(i,j)-u(i+1,j)
                Qx_p(i,j,2)=u(i+1,j)+(w*(b0-b1)+b1)/2.d0
                !
                b0=v(i+1,j)-v(i+2,j)
                b1=v(i,j)-v(i+1,j)
                Qx_p(i,j,3)=v(i+1,j)+(w*(b0-b1)+b1)/2.d0
                !
                b0=p(i+1,j)-p(i+2,j)
                b1=p(i,j)-p(i+1,j)
                Qx_p(i,j,4)=p(i+1,j)+(w*(b0-b1)+b1)/2.d0
            enddo
        enddo
    case(13,14) !
        g1=1.d0/4.d0
        g2=13.d0/12.d0
        do j=1,gridsy
            do i=0,gridsx
                bm1=share_f(i-1,j)-share_f(i-2,j)
                b0=share_f(i,j)-share_f(i-1,j)
                b1=share_f(i+1,j)-share_f(i,j)
                b2=share_f(i+2,j)-share_f(i+1,j)
            
                f10=3.d0*b0-bm1
                f11=b0+b1
                f12=3.d0*b1-b2
                f20=b0-bm1
                f21=b1-b0
                f22=b2-b1

                beta0=g1*f10*f10+g2*f20*f20+eps
                beta1=g1*f11*f11+g2*f21*f21+eps
                beta2=g1*f12*f12+g2*f22*f22+eps

                w0=1.d0/beta0
                w1=6.d0/beta1
                w=w0+w1+3.d0/beta2
                w0=w0/w
                w1=w1/w
                !                
                bm1=d(i-1,j)-d(i-2,j)
                b0=d(i,j)-d(i-1,j)
                b1=d(i+1,j)-d(i,j)
                b2=d(i+2,j)-d(i+1,j)
            
                q0=d(i,j)-bm1/3.d0+5.d0*b0/6.d0
                q1=d(i,j)+b1/3.d0+b0/6.d0
                q2=d(i,j)+2.d0*b1/3.d0-b2/6.d0
                
                Qx_m(i,j,1)=q2+w0*(q0-q2)+w1*(q1-q2)
                !
                bm1=u(i-1,j)-u(i-2,j)
                b0=u(i,j)-u(i-1,j)
                b1=u(i+1,j)-u(i,j)
                b2=u(i+2,j)-u(i+1,j)
            
                q0=u(i,j)-bm1/3.d0+5.d0*b0/6.d0
                q1=u(i,j)+b1/3.d0+b0/6.d0
                q2=u(i,j)+2.d0*b1/3.d0-b2/6.d0
                
                Qx_m(i,j,2)=q2+w0*(q0-q2)+w1*(q1-q2)
                !
                bm1=v(i-1,j)-v(i-2,j)
                b0=v(i,j)-v(i-1,j)
                b1=v(i+1,j)-v(i,j)
                b2=v(i+2,j)-v(i+1,j)
            
                q0=v(i,j)-bm1/3.d0+5.d0*b0/6.d0
                q1=v(i,j)+b1/3.d0+b0/6.d0
                q2=v(i,j)+2.d0*b1/3.d0-b2/6.d0
                
                Qx_m(i,j,3)=q2+w0*(q0-q2)+w1*(q1-q2)
                !
                bm1=p(i-1,j)-p(i-2,j)
                b0=p(i,j)-p(i-1,j)
                b1=p(i+1,j)-p(i,j)
                b2=p(i+2,j)-p(i+1,j)
            
                q0=p(i,j)-bm1/3.d0+5.d0*b0/6.d0
                q1=p(i,j)+b1/3.d0+b0/6.d0
                q2=p(i,j)+2.d0*b1/3.d0-b2/6.d0
                
                Qx_m(i,j,4)=q2+w0*(q0-q2)+w1*(q1-q2)
    
                bm1=share_f(i+2,j)-share_f(i+3,j)
                b0=share_f(i+1,j)-share_f(i+2,j)
                b1=share_f(i,j)-share_f(i+1,j)
                b2=share_f(i-1,j)-share_f(i,j)
            
                f10=3.d0*b0-bm1
                f11=b0+b1
                f12=3.d0*b1-b2
                f20=b0-bm1
                f21=b1-b0
                f22=b2-b1

                beta0=g1*f10*f10+g2*f20*f20+eps
                beta1=g1*f11*f11+g2*f21*f21+eps
                beta2=g1*f12*f12+g2*f22*f22+eps

                w0=1.d0/beta0
                w1=6.d0/beta1
                w=w0+w1+3.d0/beta2
                w0=w0/w
                w1=w1/w
                !
                bm1=d(i+2,j)-d(i+3,j)
                b0=d(i+1,j)-d(i+2,j)
                b1=d(i,j)-d(i+1,j)
                b2=d(i-1,j)-d(i,j)
            
                q0=d(i+1,j)-bm1/3.d0+5.d0*b0/6.d0
                q1=d(i+1,j)+b1/3.d0+b0/6.d0
                q2=d(i+1,j)+2.d0*b1/3.d0-b2/6.d0
                
                Qx_p(i,j,1)=q2+w0*(q0-q2)+w1*(q1-q2)
                !
                bm1=u(i+2,j)-u(i+3,j)
                b0=u(i+1,j)-u(i+2,j)
                b1=u(i,j)-u(i+1,j)
                b2=u(i-1,j)-u(i,j)
            
                q0=u(i+1,j)-bm1/3.d0+5.d0*b0/6.d0
                q1=u(i+1,j)+b1/3.d0+b0/6.d0
                q2=u(i+1,j)+2.d0*b1/3.d0-b2/6.d0
                
                Qx_p(i,j,2)=q2+w0*(q0-q2)+w1*(q1-q2)
                !
                bm1=v(i+2,j)-v(i+3,j)
                b0=v(i+1,j)-v(i+2,j)
                b1=v(i,j)-v(i+1,j)
                b2=v(i-1,j)-v(i,j)
            
                q0=v(i+1,j)-bm1/3.d0+5.d0*b0/6.d0
                q1=v(i+1,j)+b1/3.d0+b0/6.d0
                q2=v(i+1,j)+2.d0*b1/3.d0-b2/6.d0
                
                Qx_p(i,j,3)=q2+w0*(q0-q2)+w1*(q1-q2)
                !
                bm1=p(i+2,j)-p(i+3,j)
                b0=p(i+1,j)-p(i+2,j)
                b1=p(i,j)-p(i+1,j)
                b2=p(i-1,j)-p(i,j)
            
                q0=p(i+1,j)-bm1/3.d0+5.d0*b0/6.d0
                q1=p(i+1,j)+b1/3.d0+b0/6.d0
                q2=p(i+1,j)+2.d0*b1/3.d0-b2/6.d0
                
                Qx_p(i,j,4)=q2+w0*(q0-q2)+w1*(q1-q2)
            enddo
        enddo 
    case(16) !weno3is+critical point q=2 +compact third-order
        do j=1,gridsy
            i=0
            b0=share_f(i,j)-share_f(i-1,j)
            b1=share_f(i+1,j)-share_f(i,j)
            b2=share_f(i+2,j)-share_f(i+1,j)
            is2=4.d0*(b1-b0)**2
            beta2=(dabs(b1+b2)-dabs(3.d0*b2-b1))**2
            theta=dabs(is2-beta2)
            if(theta>min(is2,beta2)) then
                beta0=(dabs(b0+b1)-dabs(3.d0*b0-b1))**2+eps
                beta1=(dabs(b0+b1)-dabs(3.d0*b1-b0))**2+eps
                w=1.d0/(1.d0+2.d0*(beta0/beta1)**2)
                !
                b0=d(i,j)-d(i-1,j)
                b1=d(i+1,j)-d(i,j)
                Qx_m(i,j,1)=d(i,j)+(w*(b0-b1)+b1)/2.d0
                !
                b0=u(i,j)-u(i-1,j)
                b1=u(i+1,j)-u(i,j)
                Qx_m(i,j,2)=u(i,j)+(w*(b0-b1)+b1)/2.d0
                !
                b0=v(i,j)-v(i-1,j)
                b1=v(i+1,j)-v(i,j)
                Qx_m(i,j,3)=v(i,j)+(w*(b0-b1)+b1)/2.d0
                !
                b0=p(i,j)-p(i-1,j)
                b1=p(i+1,j)-p(i,j)
                Qx_m(i,j,4)=p(i,j)+(w*(b0-b1)+b1)/2.d0
            else
                b0=d(i,j)-d(i-1,j)
                b1=d(i+1,j)-d(i,j)
                Qx_m(i,j,1)= d(i,j)+(2.d0*b1+b0)/6.d0
                b0=u(i,j)-u(i-1,j)
                b1=u(i+1,j)-u(i,j)
                Qx_m(i,j,2)= u(i,j)+(2.d0*b1+b0)/6.d0
                b0=v(i,j)-v(i-1,j)
                b1=v(i+1,j)-v(i,j)
                Qx_m(i,j,3)= v(i,j)+(2.d0*b1+b0)/6.d0
                b0=p(i,j)-p(i-1,j)
                b1=p(i+1,j)-p(i,j)
                Qx_m(i,j,4)= p(i,j)+(2.d0*b1+b0)/6.d0
            endif
            do i=1,gridsx
                b0=share_f(i,j)-share_f(i-1,j)
                b1=share_f(i+1,j)-share_f(i,j)
                b2=share_f(i+2,j)-share_f(i+1,j)
                is2=4.d0*(b1-b0)**2
                beta2=(dabs(b1+b2)-dabs(3.d0*b2-b1))**2
                theta=dabs(is2-beta2)
                beta0=(dabs(b0+b1)-dabs(3.d0*b0-b1))**2+eps
                beta1=(dabs(b0+b1)-dabs(3.d0*b1-b0))**2+eps
                if( theta>min(is2,beta2).and.dabs(beta0-beta1)>min(beta0,beta1) ) then
                    w=1.d0/(1.d0+2.d0*(beta0/beta1)**2)
                    !
                    b0=d(i,j)-d(i-1,j)
                    b1=d(i+1,j)-d(i,j)
                    Qx_m(i,j,1)=d(i,j)+(w*(b0-b1)+b1)/2.d0
                    !
                    b0=u(i,j)-u(i-1,j)
                    b1=u(i+1,j)-u(i,j)
                    Qx_m(i,j,2)=u(i,j)+(w*(b0-b1)+b1)/2.d0
                    !
                    b0=v(i,j)-v(i-1,j)
                    b1=v(i+1,j)-v(i,j)
                    Qx_m(i,j,3)=v(i,j)+(w*(b0-b1)+b1)/2.d0
                    !
                    b0=p(i,j)-p(i-1,j)
                    b1=p(i+1,j)-p(i,j)
                    Qx_m(i,j,4)=p(i,j)+(w*(b0-b1)+b1)/2.d0
                else
                    b1=d(i+1,j)-d(i,j)
                    Qx_m(i,j,1)= 1.5d0*d(i,j)+b1/4.d0-Qx_m(i-1,j,1)/2.d0
                    b1=u(i+1,j)-u(i,j)
                    Qx_m(i,j,2)= 1.5d0*u(i,j)+b1/4.d0-Qx_m(i-1,j,2)/2.d0
                    b1=v(i+1,j)-v(i,j)
                    Qx_m(i,j,3)= 1.5d0*v(i,j)+b1/4.d0-Qx_m(i-1,j,3)/2.d0
                    b1=p(i+1,j)-p(i,j)
                    Qx_m(i,j,4)= 1.5d0*p(i,j)+b1/4.d0-Qx_m(i-1,j,4)/2.d0
                endif
            enddo     
            i=gridsx
            b0=share_f(i+1,j)-share_f(i+2,j)
            b1=share_f(i,j)-share_f(i+1,j)
            b2=share_f(i-1,j)-share_f(i,j)
            is2=4.d0*(b1-b0)**2
            beta2=(dabs(b1+b2)-dabs(3.d0*b2-b1))**2
            theta=dabs(is2-beta2)
            if(theta>min(is2,beta2)) then
                beta0=(dabs(b0+b1)-dabs(3.d0*b0-b1))**2+eps
                beta1=(dabs(b0+b1)-dabs(3.d0*b1-b0))**2+eps
                w=1.d0/(1.d0+2.d0*(beta0/beta1)**2)
                !
                b0=d(i+1,j)-d(i+2,j)
                b1=d(i,j)-d(i+1,j)
                Qx_p(i,j,1)=d(i+1,j)+(w*(b0-b1)+b1)/2.d0
                !
                b0=u(i+1,j)-u(i+2,j)
                b1=u(i,j)-u(i+1,j)
                Qx_p(i,j,2)=u(i+1,j)+(w*(b0-b1)+b1)/2.d0
                !
                b0=v(i+1,j)-v(i+2,j)
                b1=v(i,j)-v(i+1,j)
                Qx_p(i,j,3)=v(i+1,j)+(w*(b0-b1)+b1)/2.d0
                !
                b0=p(i+1,j)-p(i+2,j)
                b1=p(i,j)-p(i+1,j)
                Qx_p(i,j,4)=p(i+1,j)+(w*(b0-b1)+b1)/2.d0
            else
                b0=d(i+1,j)-d(i+2,j)
                b1=d(i,j)-d(i+1,j)
                Qx_p(i,j,1)=d(i+1,j)+(2.d0*b1+b0)/6.d0    
                !
                b0=u(i+1,j)-u(i+2,j)
                b1=u(i,j)-u(i+1,j)
                Qx_p(i,j,2)=u(i+1,j)+(2.d0*b1+b0)/6.d0    
                !
                b0=v(i+1,j)-v(i+2,j)
                b1=v(i,j)-v(i+1,j)
                Qx_p(i,j,3)=v(i+1,j)+(2.d0*b1+b0)/6.d0    
                !
                b0=p(i+1,j)-p(i+2,j)
                b1=p(i,j)-p(i+1,j)
                Qx_p(i,j,4)=p(i+1,j)+(2.d0*b1+b0)/6.d0    
            endif  
            do i=gridsx-1,0,-1
                b0=share_f(i+1,j)-share_f(i+2,j)
                b1=share_f(i,j)-share_f(i+1,j)
                b2=share_f(i-1,j)-share_f(i,j)  
                is2=4.d0*(b1-b0)**2
                beta2=(dabs(b1+b2)-dabs(3.d0*b2-b1))**2
                theta=dabs(is2-beta2)
                beta0=(dabs(b0+b1)-dabs(3.d0*b0-b1))**2+eps
                beta1=(dabs(b0+b1)-dabs(3.d0*b1-b0))**2+eps
                if( theta>min(is2,beta2).and.dabs(beta0-beta1)>min(beta0,beta1) ) then
                    w=1.d0/(1.d0+2.d0*(beta0/beta1)**2)
                    !
                    b0=d(i+1,j)-d(i+2,j)
                    b1=d(i,j)-d(i+1,j)
                    Qx_p(i,j,1)=d(i+1,j)+(w*(b0-b1)+b1)/2.d0
                    !
                    b0=u(i+1,j)-u(i+2,j)
                    b1=u(i,j)-u(i+1,j)
                    Qx_p(i,j,2)=u(i+1,j)+(w*(b0-b1)+b1)/2.d0
                    !
                    b0=v(i+1,j)-v(i+2,j)
                    b1=v(i,j)-v(i+1,j)
                    Qx_p(i,j,3)=v(i+1,j)+(w*(b0-b1)+b1)/2.d0
                    !
                    b0=p(i+1,j)-p(i+2,j)
                    b1=p(i,j)-p(i+1,j)
                    Qx_p(i,j,4)=p(i+1,j)+(w*(b0-b1)+b1)/2.d0
                else
                    b1=d(i,j)-d(i+1,j)
                    Qx_p(i,j,1)=1.5d0*d(i+1,j)+b1/4.d0-Qx_p(i+1,j,1)/2.d0 
                    !
                    b1=u(i,j)-u(i+1,j)
                    Qx_p(i,j,2)=1.5d0*u(i+1,j)+b1/4.d0-Qx_p(i+1,j,2)/2.d0  
                    !
                    b1=v(i,j)-v(i+1,j)
                    Qx_p(i,j,3)=1.5d0*v(i+1,j)+b1/4.d0-Qx_p(i+1,j,3)/2.d0  
                    !
                    b1=p(i,j)-p(i+1,j)
                    Qx_p(i,j,4)=1.5d0*p(i+1,j)+b1/4.d0-Qx_p(i+1,j,4)/2.d0 
                endif
            enddo 
        enddo
    case default
        stop'share_space_x: zero pivot at first row'
    end select
    
end subroutine share_space_x